something on the similar lines(in terms of enhancing tournament algorithms)-- we can also consider taking this reusing the pool of workers approach maybe for tournament_is_strongly_connected (in line 80) -- and that might show some more improvements in terms of speedups-- but all our functions are currently wrapped in a _configure_if_nx_active-- which is essentially a joblib.parallel_config context-- so this would require us to change how _configure_if_nx_active works currently.

PS: this doesn't have to be handled in this PR-- just something we can think about in the long run.

Hi @Schefflera-Arboricola and @dschult, I'm trying to run asv benchmarks for is_reachable from my branch, but I keep running into a numpy import error. Can you help me out with this? I suspect it might be because NumPy isn't a default dependency. I tried adding it to pyproject.toml, but that didn't resolve the issue.

How is your environment set up? "conda/mamba" or "venv + pip"?
What script/incantation are you using? Can you give me enough info to recreate the issue locally? (I have your branch locally)

Apologies for the oversight! I was able to figure it out from the documentation. Thank you @dschult!

I tried out the following benchmarks:

1. Compared networkx's recent change to the numpy implementation of is_reachable() on this branch.
   Results:

========== ============ ========== =========== ============ ============
--                                   num_nodes                          
---------- -------------------------------------------------------------
backend        50         100         200         400          800     
========== ============ ========== =========== ============ ============
parallel   6.19±0.4ms   14.5±2ms    44.7±1ms    166±1ms      655±2ms   
networkx    18.5±3ms    74.3±3ms   288±0.7ms   1.13±0.01s   4.53±0.03s 
========== ============ ========== =========== ============ ============

2. Compared networkx's recent change to a parallel version of the same (i.e, parallelised the pure Python Networkx implementation without Numpy)
   Results:

========== ============ ========== ============ ========== ============
--                                  num_nodes                          
---------- ------------------------------------------------------------
backend        50         100         200         400         800     
========== ============ ========== ============ ========== ============
parallel   1.15±0.1ms   5.00±1ms   17.5±0.3ms   67.8±1ms   270±0.7ms  
networkx    18.3±1ms    73.7±1ms    289±3ms     1.12±0s    4.59±0.03s 
========== ============ ========== ============ ========== ============

The pure Python implementation performs better. I haven't committed the pure Python code yet because I wanted to confirm a few things first. Do these benchmarks imply that the process of converting a graph into a numpy matrix is more expensive than that of passing the huge graph to each core? Because n_jobs is None by default, for these benchmarks, I'm unsure how we are seeing the above speedups.

Pure Python parallelism after the recent Networkx commit, for reference:

def two_neighborhood_close(G, chunk):
        for v in chunk:
            v_adj = G._adj[v]
            S = {
                x
                for x, x_pred in G._pred.items()
                if x == v or x in v_adj or any(z in v_adj for z in x_pred)
            }
            if s in S and t not in S and is_closed(G, S):
                return True
        return False

    def is_closed(G, S):
        return all(u in S or all(v in unbrs for v in S) for u, unbrs in G._adj.items())

    if hasattr(G, "graph_object"):
        G = G.graph_object

    n_jobs = nxp.get_n_jobs()

    if get_chunks == "chunks":
        node_chunks = nxp.chunks(G, n_jobs)
    else:
        node_chunks = get_chunks(G)

    return not any(
        Parallel()(delayed(two_neighborhood_close)(G, chunk) for chunk in node_chunks)
    )

Do these benchmarks imply that the process of converting a graph into a numpy matrix is more expensive than that of passing the huge graph to each core? Because n_jobs is None by default, for these benchmarks, I'm unsure how we are seeing the above speedups.

Which version of config is used when n_jobs is None? That is, how many jobs are you actually running? The last benchmark with parallel w/ python vs sequential python looks like parallel is about 16 times faster. So that is either 16 cores or 4 cores with quadratic speed-up. I don't see any quadratic speedup looking at the code. But it is possible I am missing something.

Can you point me to your benchmarking code? Is it already merged in nx-parallel, or is it only on your local machine?

That is, how many jobs are you actually running?

I didn’t modify the number of jobs anywhere in the codebase, so by default, n_jobs=None is being used by joblib which would mean 1 core. I also verified this by checking the activity monitor while running the benchmark — the CPU usage confirmed that only a single core was active during the run.

Can you point me to your benchmarking code? Is it already merged in nx-parallel, or is it only on your local machine?

The local changes I made were:

• replacing

   _ = nx.tournament.is_reachable(G, 1, num_nodes, backend=backend)
  

  with:

  _ = nx.tournament.is_reachable(G, 0, num_nodes - 1, backend=backend)
  

  to prevent the target node from going out of bounds.
• adding "networkx" to the list of backends in common.py.

I ran the timing script on both numpy implementation and the pure python implementation.
Numpy:

200
0.16524791717529297
0.0043315887451171875
Finished <function is_reachable at 0x1037ba980>
400
0.04875922203063965
0.017219066619873047
Finished <function is_reachable at 0x1037ba980>
800
0.1309349536895752
0.06830096244812012
Finished <function is_reachable at 0x1037ba980>
1600
0.4890291690826416
0.31351304054260254
Finished <function is_reachable at 0x1037ba980>
is_reachable

Python:

200
0.2303180694580078
0.004209995269775391
Finished <function is_reachable at 0x105c5a980>
400
0.056372880935668945
0.016617774963378906
Finished <function is_reachable at 0x105c5a980>
800
0.30538177490234375
0.06642794609069824
Finished <function is_reachable at 0x105c5a980>
1600
1.9619712829589844
0.3054389953613281
Finished <function is_reachable at 0x105c5a980>
is_reachable

These timings clearly show that the NumPy version outperforms the pure Python parallel version as the graph size grows. This is mainly because converting the graph to an adjacency matrix is a one-time cost, while passing the entire Graph object across multiple cores in the pure Python version repeatedly triggers expensive pickling and unpickling. The current implementation avoids this by storing the matrix on disk and allowing each worker to read it directly (using mem-mapping). I’m happy with this approach instead of converting it to pure Python parallel code, but please let me know if you have any other suggestions.

I believe this PR is ready if there aren’t any other specifics to address :)